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When auxiliary information is available at the design stage, samples may be selected by means of 
balanced sampling. Deville and Tille proposed in 2004 a general algorithm to perform balanced 
sampling, named the cube method. In this paper, we are interested in a particular case of the cube 
method named pivotal sampling, and first described by Deville and Tille in 1998. We show that 
this sampling algorithm, when applied to units ranked in a fixed order, is equivalent to Deville's 
systematic sampling, in the sense that both algorithms lead to the same sampling design. This 
characterization enables the computation of the second-order inclusion probabilities for pivotal 
sampling. We show that the pivotal sampling enables to take account of an appropriate ordering 
of the units to achieve a variance reduction, while limiting the loss of efficiency if the ordering 
is not appropriate. 

Keywords: balanced sampling; cube method; design effect; sampling algorithm; second order 
inclusion probabilities; unequal probabilities 

1. Introduction 

When auxiliary information is available at the design stage, samples may be selected by 
means of balanced sampling. The variance of the Horvitz-Thompson (HT) estimator is 
then reduced, since it is approximately given by that of the residuals of the variable of 
interest on the balancing variables. Deville and Tille [6] proposed a general algorithm for 
balanced sampling, named the cube method. This sampling algorithm enables the selection 
of balanced samples with any number of balancing variables, and any prescribed set of 
inclusion probabilities. 

In order to measure the gain in efficiency provided by the cube method, Deville and 
Tille [7] proposed several variance approximations. They suppose that the sampling de- 
sign is exactly balanced, and performed with maximum entropy among sampling designs 
balanced on the same balancing variables, with the same inclusion probabilities. Then, 
under an additional assumption of asymptotic normality of the multivariate HT-cstimator 
under Poisson sampling, the variance approximations are derived. The assumption of ex- 
act balancing may be closely respected, if the number of balancing variables remains 
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small with regard to the sample size; otherwise, the balancing error must be taken into 
account in variance estimation, see Breidt and Chauvet [9]. The second assumption is 
related to the entropy of the sampling design: the variance approximations proposed by 
Deville and Tille [7] are unlikely to hold if this assumption is not satisfied. 

A practical way to increase the entropy of a sampling design is to sort the popu- 
lation randomly before sampling. However, this preliminary randomization step is not 
systematically included in the sampling process. This is a common practice to sort the 
population with respect to some auxiliary variable before the sampling, so as to ben- 
efit from a stratification effect. In France, Census surveys are conducted annually; the 
detailed methodology is described in Godinot [8]. Each large municipality (10 000 inhab- 
itants or more in 1999) is the subject of an independent sampling design and is stratified 
according to the type of address (large addresses, new addresses, or other addresses). 
In each stratum, the addresses are divided into 5 rotation groups. Each year, all the 
addresses within one rotation group (for the strata of large addresses and new addresses) 
or within a sub-sample (for the stratum of other addresses) are surveyed. In the stratum 
of other addresses, the sub-sample is obtained by first, sorting the addresses with respect 
to the descending number of dwellings, and then, applying the cube method. In such 
cases, the conditions for the variance approximations proposed by Deville and Tille [7] 
to hold are clearly not respected. 

We are interested in a particular case of the cube method, called pivotal sampling 
(Deville and Tille [5]), obtained when the only balancing condition is given by the variable 
of inclusion probabilities. That is, the cube method with the sole fixed-size constraint 
amounts to pivotal sampling. This algorithm is an exact sampling procedure, which 
respects a prescribed set of inclusion probabilities, is strictly without replacement and 
leads to fixed-size designs. In this paper, we show that the pivotal sampling algorithm, 
when applied to units ranked in a fixed order, is equivalent to an algorithm proposed in 
Deville [4], and known in the literature as Deville 's systematic sampling (Tille [13]). The 
two algorithms are equivalent, in the sense that both lead to the same sampling design. 
In particular, the computation of the second-order inclusion probabilities developed in 
Deville [4] may be readily applied to pivotal sampling. This provides an answer to a 
problem raised by Bondesson and Grafstrom [1], page 7. Devillc's systematic sampling 
has similarities with Markov chain designs introduced by Breidt [2]. It has found uses in 
the context of longitudinal surveys, see Nedyalkova, Qualite, and Tille [11]. 

The paper is organized as follows. In Section 2, the notation is defined. Ordered pivotal 
sampling and Deville's systematic sampling are presented in Sections 3 and 4, respec- 
tively, and some useful results are derived. The second-order inclusion probabilities for 
ordered pivotal sampling are given in Section 5. Some results which illustrate the practical 
interest of ordered pivotal sampling arc presented in Section 6. 

2. Notation 

Consider a finite population U consisting of N sampling units that may be represented 
by integers k = 1, . . . , N . We assume that the order of the units in the population is 
fixed prior to sampling, and may be confounded with the natural order of their in- 
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Figure 1. Inclusion probabilities and cross-border units in microstratum Ui, for population U. 



dexes. A sample s, denned as a subset of U, is selected with inclusion probabilities 
7r = (tti , . . . , 7Tjv)'- We assume without loss of generality that < Hk < 1 for any unit k in 
U, with n = X^fcec/ ^ ^ ne sam pl e size. Let tt^i denote the probability that units fc and I 
are selected jointly in the sample. 

We define Vk = X^zLi n l f° r an y un ^ k & U , with Vb = 0. A unit k is said to be cross- 
border if Vfc_i < i and Vk > i for some nonnegative integer i. The cross-border units 
are denoted as fcj,i = l,...,n — 1, and we note Oj = i — Vj^-i and fr; = Vfe; — i. The 
microstratum C/j, i = 1, . . . , n, is defined as 

?7i = {k e ?7; fci-i < < fci}, (2.1) 

with fco = and k n = N+ 1. To fix ideas, useful quantities for population (7 are presented 
in Figure 1. 

The microstrata are generally overlapping, since one cross-border unit may belong to 
two adjacent microstrata: the cross-border unit ki belongs both to the microstratum Ui 
(with an associated probability a{) and to the microstratum C/j+i (with an associated 
probability bi). In the particular case when = i, we have bi = 0. To avoid the introduc- 
tion of specific notations for such cases, we consider in Sections 2-5 that, in this situation, 
the cross-border unit ki belongs to the microstratum U^i as a "phantom unit," that 
is, with an associated probability equal to 0. In Section 6, we simply consider that the 
cross-border unit ki belongs to the microstratum Ui only in such situations. 

The N sampling units are grouped to obtain a population U c = {ui, . . . , U2 n -i} of 
clusters. There are the clusters of cross-border units (n — 1 singletons), denoted as un 
with associated probability (f>2i = i^ki for i = 1, . . . ,n — 1. There are the n clusters of 
units that are not cross-borders and that are between two consecutive integers, denoted 
as ii2i-\ with associated probability (p2i-i = Vki-i ~ ^4i-u f° r * — 1, We note 

ijj = (01, . . . , 4>2n-i)' ■ To fix ideas, useful quantities for population U c are presented in 
Figure 2. If (at least) one of the cross-border units in Ui has a large inclusion probability, 
there may not exist any non-cross-border unit between integers i — 1 and i, so that the 
cluster U2i-i is empty. To avoid the need for specific notations for such cases, we may 
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Figure 2. Inclusion probabilities and cross-border units in microstrata Ui and Ui+i for popu- 
lation U c . 

view this situation as a particular case of our framework by allowing a cluster u%i— i to 
be a "phantom cluster," that is, an empty cluster with associated probability <f>2i-i equal 
to 0. For example, suppose that iV = 8, n = 4 and tt = (0.2, 0.5, 0.3, 0.4, 0.9, 0.8, 0.5, 0.4)'. 
We obtain the 4 microstrata U x = {1,2,3}, U 2 = {3,4,5}, U 3 = {5,6} and C/ 4 = {6,7,8}. 
In particular, we have Oi = 0.3 = and b\ = 0, so that the cross-border unit 3 is a 
phantom unit for the microstratum C/jj. Also, we obtain 7 clusters (see Table 1): the 
cluster M5 is empty, with an associated probability equal to zero. 

3. Ordered pivotal sampling 

A general algorithm for pivotal sampling is described in Deville and Tille [5]. In the 
version presented in Algorithm 1, the order of the sampling units is explicitly taken into 
account. We call it ordered pivotal sampling to avoid confusion. At each step, one or two 
coordinates of ir(t) are randomly rounded to or 1, and remain there forever. In at most 
N steps, the final sample is obtained. 

Roughly speaking, the algorithm may be summarized as follows. At the beginning, 
in microstratum U\ (i = 1), the two first units 1 and 2 fight, the loser is definitely 
eliminated while the survivor (denoted as Jo) gets the sum of their probabilities and 
then faces the following unit. The fights go on until the accumulated probability exceeds 
1 , which occurs at time t = ki when the survivor Jq faces the cross-border unit ki ■ One 
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Algorithm 1 Ordered Pivotal Sampling with parameter n 

1. We initialize with i = 1, Jo = 1 and tt(0) = 7T. 

2. For t = 2,...,N, do: 

(a) If ?n G C7 \ { J;_i,i}, then 7r rn (/j) = 7r TO (< — 1). 

(b) If 7rj i _ 1 (t - 1) + * t (t- 1) < 1, let Ai(t) = ^-^Effe^j . Then 

i. with probability Ai(t), let 

[nj t -dt),ir t (fii\ = [TJ*_i(*-l) + Tt(t - 1),0]; 

ii. with probability 1 — Ai(i), let J^_i = t and 

K/i-i (*),**(*)] = [O.TTj^t - 1) +7T t (t- 1)]. 

(c) If 7Tj i _ 1 (t - 1) + n(t - 1) > 1, let Ai(t) = a-^^'ffxV-^t-D - Then 

i. with probability Ai(t), let Wi = Ji-i, let Jj =i and 

[TTj^^tJ.TTt^)] - [l.TTj^t - 1) +7T t (t - 1) - 1]; 

ii. with probability 1 — Ai(t), let IVj = t, let = Jj_i and 

[Trj^^tJ.TrtW] - [^^(t - 1) + 7r t (t - 1) - 1, 1]; 

iii. let i = i + l. 

3. The sample is given by {W\, ■ ■ ■ , W n }. 



of the two remaining units, denoted as W\ , wins and is then definitely selected in the 
sample while the other one, denoted as J±, jumps to the microstratum U2- 

More generally, in microstratum Ui, the first unit fcj_i is replaced with the unit Jj_i 
which jumps from the microstratum Ui-\. The units Ji-i and fcj_i + 1 fight, the survivor 
gets the sum of their probabilities and then faces the next unit. The fights go on until 
the survivor Jj_i faces the cross-border unit ki. One of the two remaining units (Wi) 
wins and is then definitely selected in the sample while the other one (J,) jumps to the 
following microstratum. Lemma 3.1 states that Algorithm 1 may alternatively be seen 
as a two-stage procedure. The proof follows from definition, and is thus omitted. 

Lemma 3.1. Ordered pivotal sampling with parameter tt may be obtained by two-stage 
sampling, where a sample s c of n clusters is first selected in U c by means of ordered 
pivotal sampling with parameter tp, and one unit k is then selected in each Uj G s c with 
a probability proportional to 7Tfc . 

We assume that a sample S op is selected in U c by means of ordered pivotal sam- 
pling with parameter -0, and we let X\ < • • • < X n denote the units selected in the 
sample, ranked in ascending order. Lemma 3.2 states useful relations between on the one 
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hand, the sampled units Xi, and on the other hand, the winners Wi and jumpers Ji. 
Lemma 3.3 gives the probabilities for the different outcomes in the case of a non-cross- 
border unit M2-8-1- 

Lemma 3.2. In case of ordered pivotal sampling with parameter ip, we have 

{X, = u 2l _ 2 } => { J,_! G {X u . . . , X t }}, (3.1) 
{Xi = uai-i} =► {W t = u 2i -i} U {Ji = ujh-i}, (3.2) 
{X t = u 2l } => {J, i {X 1 , . . .,Xi}}. (3.3) 

Proof. Assume that Xi = u 2 j_ 2 . This implies that i units exactly are selected in the 
i — 1 first microstrata t/i, . . . , f/j-i. On the other hand, if Jj_i ^ {-X"i, . . . , ATi} the unit 
Jj_i is not selected in the sample so that at most i — 1 units are selected in U±, . . . , J7f-i. 
This proves (3.1), and by a similar argument we obtain (3.3). It is easily seen that (3.2) 
holds, since the selection of u 2 i— l implies that this unit is either the winner Wi or the 
jumper Ji in the microstratum L/j. □ 

Lemma 3.3. In case of ordered pivotal sampling with parameter ip, we have 

(rxr s (l-Oi-6i_i)(l-Oi-6i) 
pr(Wi = u 2i _i = — — , 3.4 

(1 - Oj)(l - Oj) 

/ t \ Oi(l-Oi-6i_i) , 
pr(J t = u 2l _i = — — — , (3.5 

(1 - Oi)(l - &i) 

pr(Xi = uat_i) = 1 — Qi — 6i_i. (3.6) 

Proof. The event 

{W l ^u 2l -i} 

may be alternatively interpreted as follows: in the fight between Jj_i and u 2 i-i, the unit 
u%i-x survives; then in the next fight, the unit it 2 j_i is the selected unit Wi, while the 
unit un is the jumping unit Jj. Consequently, we have: 

pr(Wj = tt 2 *-i) = : x — — - — , 

1 - di l — Oi 

which gives (3.4). Similarly, we obtain 

, T s l-6 ? -i -a.; a-i 

pr{J. t = u 2i -i) = x -, 

1 - a, 1 - bi 

which gives (3.5). We now consider equation (3.6). Since 



{X t = ti 2l -i} => G 5 op } 
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and 

pr(u 2i _i e S op ) = 1 - a, - bi_i, 

it suffices to show that 

{u 2 i-i eS op }=>{Xi = U2i-i}- (3.7) 

Since {u 2 i-i € S^p} implies that m 2 ;-i survives in its duel against Ji_i, this in turn 
implies that Jj_i ^ {-Xi, . . . , X^}. In other words, \u<&-\ S Sop} implies that exactly i — 1 
units smaller than un-\ were selected, which proves (3.7). □ 

Finally, let U Cii = {uai-2, • • • , «2n-i}> V>i = 02»-i> • ••)^i>- • -,</>2n-i)', and S oP) j be 

a random sample selected in U c .i by means of ordered pivotal sampling with parameter 
ipi. Lemma 3.4 establishes some relations for conditional inclusion probabilities in S oPi i 
of the first units in U c j. 

Lemma 3.4. 

pr{u 2l £ S op ,i,U2i-i i- S p^\u 2 i~2 & S p }l ) 
b; 



(3.8) 



1 — Cti ' 

pr(u 2 i+i 6 S op ,i,u 2 i £ S op ,i,U2i-i i- S op ,i\u2i-2 e S op j) 

_ (1 - Qj - - &j - Oi+l) 

(l-Oi)(l-6i) 

u 2 i-i S oP}l \u 2i - 2 6 Sop,z) 

_ (1 - Oj - bpOi-j-i 
(1-00(1-60 ■ 



(3.9) 



(3.10) 



Proof. To fix ideas, the first units in population U c ^ and related quantities are presented 
in Figure 3. 

We first consider equation (3.8). Since b^± is the first-order inclusion probability of 
unit «2i-2 in sample S op ^, we have 

pr(u2i-2 €S oPti ) = 6j_i. (3-11) 

On the other hand, the event 

{u 2 i € S op ^,u 2 i^i f. S op .i,u 2 i- 2 £ S oPt i} 

may be alternatively interpreted as follows: in the first fight, the unit u 2 i- 2 survives 
against the unit u 2 i-ii in the second fight, any of the two units u 2 i- 2 or u 2 i is the 
selected unit W\, while the other is the jumping unit J\\ then, the jumping unit J\ is 
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Figure 3. Inclusion probabilities and cross-border units in the two first microstrata of popula- 
tion U c , i ■ 



selected during one of the following fights. Consequently, we have: 
pr(u 2i € S op .i,u 2l -i £ S op ,i,u 2 i-2 € 



(3.12) 

x 1 x bi, 



1 - Oi 



and equation (3.8) follows from (3.11) and (3.12). We now consider equation (3.9). The 
event 

{u 2 i+l £ S op j, U 2 i £ S op .i, U 2 i-\ f. S OPl i, U 2 i- 2 £ S op .i} 

may be interpreted as follows: in the first fight, the unit u 2 i- 2 survives against the unit 
U2i—i', in the second fight, u 2 i- 2 is the selected unit W\, while u 2 i is the jumping unit J\\ 
in the third fight, the unit M2i+i survives against the unit u 2 i\ then, the unit 1*21+1 is 
selected during one of the following fights. Consequently, we have: 

V r ( u 2i+1 € S op ,i,U 2 i f: S OPt i,U 2 i-l f. S OP: i,U2i-2 S S oP: i) 

= - x— x— x(l-a l+1 ) 3.13) 

1 - Oj 1 - bi 1 - a i+ i 

_ 6j— i(l - a-i - bj){l - h - fflj+i) 
(l-ai)(l-bi) 

which, together with (3.11), leads to (3.9). Finally, we consider equation (3.10). The event 

{u 2 i+ 2 £ S OPt i, M2i+1 ^ S OPy i,U 2 i £ S op ,i,U 2 i-l S OPi i,U 2 i- 2 £ S op> i} 

may be interpreted as follows: in the first fight, the unit u 2 i~ 2 survives against the unit 
u 2 i-\] in the second fight, U2i-2 is the selected unit W\, while u 2 i is the jumping unit J\\ 
in the third fight, any of the two units Jj = u 2 i or U2i+i survives; in the fourth fight, 
«2i+2 is the selected unit W 2 , while the other unit is the jumper J 2 ; then, the unit J 2 is 
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not selected during one of the following fights. Consequently, we have: 



pr{u2i+2 & S OPi i,U2i+l f. S OP: i,U2i f. S OPi i,U2i-l ^ S op 4,U2i-2 S S op ,i) 




(3.14) 



which gives (3.10). 



□ 



4. Deville's systematic sampling 

The sampling algorithm known in the literature as Deville's systematic sampling (Dcv- 
illc [4]; Title [13]) is presented in Algorithm 2. This algorithm proceeds in n sub-samplings 
of size 1 in the microstrata U±, . . . , U n , and the random variables Wi which indicate the 
sampled units are generated so that a cross-border unit ki-i may not be selected twice 
in the sample: at step i, one unit denoted as Yi is drawn in Ui if was not selected 
at step i — 1, and in Ui \ {fci„i} otherwise. This sampling algorithm may be particularly 
useful in the context of business surveys, when a fine stratification is used leading to small 
and possibly non-integer sample size inside (micro)strata. Deville's systematic sampling 
directly handles the rounding problem, since any unit for which the sampling outcome 
is still undecided is moved to the next stratum, where the final sampling decision is then 
obtained. Lemma 4.1 follows from the definition of Algorithm 2. 



Algorithm 2 Deville's systematic sampling with parameter 7r 

At step 1: 

1. A distributed Uniform(0, 1) random variable w\ is generated. 

2. The unit k is selected if Vk-i <w\ < Vk- 



1. A random variable Wi is generated: 

(a) if unit ki_\ was selected at step i — 1, then Wi is generated according to a 
distributed Uniform(6j_i, 1) random variable, 

(b) otherwise, Wi is generated: 



At step i: 



• according to a distributed Uniform(0, 6,_i) random variable with proba- 
bility Oi_i6i_i{(l - Oi_i)(l - fri-i)} -1 , 

• according to a distributed Uniform(0, 1) random variable with probability 
1 -Oi-i6»-i{(l -Oi_i)(l -fej-i)} -1 . 



2. The unit k is selected if Vk-i < Wi + (i — 1) < Vk- 
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Lemma 4.1. Deville's systematic sampling with parameter tt may be obtained by two- 
stage sampling, where a sample s c ofn clusters is first selected in U c by means of Deville's 
systematic sampling with parameter ip, and one unit k is then selected in each Uj G s c 
with a probability proportional to -Kk- 

Assume that a sample is selected in U c by means of Deville's systematic sampling with 
parameter ip. The random variable i^+i which gives the result of the sampling in the 
microstratum Ui+\ only depends on the outcome of step i, so that 



The different cases for the transition probabilities in (4.1) easily follow from the definition 
of Algorithm 2, and are given below: 



P r{Y, +1 = UjlYt, ... t Yi) =pr(Y i+1 = Uj \X). 



(4.1) 



pr(Y i+1 =u J -|r 1 ,...,y;-_i,Yi =U2t-a) 



(4.2) 



1 — a,; ' 

(1 - hi - a. i+ i)(l - a l - hi) 



i = 2* + 1, 



(l-OiXl-60 
a i+ i{l - aj- bj) 

, (l-OiXl-fej) ' 



j = 2i + 2, 



pr(Y i+1 =Uj\Yi,...,Yi-i,Yi =u 2i -i) 



(4.3) 




hi) 



i = 2* + i, 



(l-Oi)(l-6 4 ) 




j = 2i + 2, 



pr(Y i+1 =tij\Yx,...,Yi-i,Yi 



U 2 i) 




(4.4) 



5. Second-order inclusion probabilities 



We can now formulate our main result. 



Theorem 5.1. Ordered pivotal sampling and Deville's systematic sampling with the 
same parameter tt induce the same sampling design. 
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Proof. From Lemmas 3.1 and 4.1, it is sufficient to prove the result in case of ordered 
systematic sampling and Deville's systematic sampling with parameter ip in the popula- 
tion U c . We only need to show that equations (4.2)-(4.4) hold in case of ordered pivotal 
sampling. Recall that we note 

U c ,i = {U2i~2, ■ ■ ■ ,U 2n -l}, 
Ipi = {bi-U<f>2i-l,-.-,<f>j,--.,(f>2n-iy, 

and that S op> i denotes a random sample selected in U Cj i by means of ordered pivotal 
sampling with parameter tpi (see Section 3). 

We first consider equation (4.2). From (3.1), we obtain: 

pr(X i+ i = u 2l \Xi,. . = u 2 i-2) 

= pr(X l+ i = u 2 i\Xi,. . .,Xi-i,Xi = U2i-2, Ji-i G {Xi, . . .,Xi}), 

which is equivalent to pr{u2i G S op ,i,U2i-i £ S oPt i\u 2 i- 2 G S oPi i), so that the result follows 
from equation (3.8). 
Similarly, we obtain 

pr(X i+ i = u 2i+ i\Xi, . . .,X l ^i,X l = u 2 i- 2 ) 

= pr(X i+ i = U2i+i\X 1 , . . ,,Xi-i,Xi = U2i-2,Ji-i G {Xi, . . .,Xi}) 

= pr(u 2 i+l G S op> i,U2i £ S op j,U 2 i-l S op ,i\ll2i-2 G S op ,i) 

_ (1 - Oj - bj)(l - bj - a i+1 ) 
(1- 00(1-60 

where the last line follows from (3.9), and 

pr(X l+ i = u 2l+2 \Xi, . . .,Xi-i,Xi = u 2 i-2) 

= pr(X i+ i = u 2i+2 \Xi, . . .,X l -i,X l = U2i-2,Ji-i G {Xi, . . .,Xi}) 

= pr{u 2 i+2 G S op ,i,U2i+l S p,i,U2i £ S op ,i,U2i-l £ S op _i\u 2 i-2 G Sop^i) 

_ (1 - dj - bi)a i+1 
~ (l-o<)(l-6i) ' 

where the last line follows from (3.10). This proves equation (4.2). The proof for equation 
(4.4) is similar, and is thus omitted. 

We now turn to equation (4.3). We introduce some further notation. Let 

U c ,i+1 = {v,2i, • ■ • , U 2n -l}, 

= (&i(02i+l)---(0ji--')02n-l) / j 
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and let S op> i+\ be a random sample selected in E/c,t+i by means of ordered pivotal sam- 
pling with parameter ipi+i- We have 

pr(X i+1 =U2i\X 1 ,...,X i _ 1 ,X i = u 2 i-i,Wi = U2i-i) 

= pr(X. l+ i =u 2 i\Xi, . ..,X i _ 1 ,X i = U2i-i,Ji = u 2l ) (5.1) 
= pr(u 2i 6 S oPii+ i) = b h 

where the second line in (5.1) comes from 

{Xi = «2i-i, Wi = u 2i -i}<=> {X t = u 2 i^i 1 J l = u 2i }. 

Also, 

pr(X i+ i =U2i\Xi,...,X i _ 1 ,X i =u 2 i-i,Ji =U2i-i) = 1, (5.2) 

since 

{Xi =u 2 i-i, J% = u 2l -i} =S> {X t = u 2i _ ll W l = u 2l } =S> {X i+ i = u 2l }. 

Further, 

pr(Wi =U2i-i\Xi,...,X i - 1 ,X i = u 2 i-i) 

= pr(Wi=U2i-i\Xi = U2i-i) (5.3) 

iy iw ■,pr(W l =u 2 i-i) 

= pr(X t =U2i-i\Wi = U2i-ij 



1 x 



pr{X r = U 2 i-l) 

(1 - oj - - a, - bj){{l - Oi)(l - h)}- 1 

I — at - 



the fourth line in (5.3) being a consequence of Lemma 3.3. The same reasoning leads to 

pr(J t =U2i-i\X 1 ,...,X i -i,X i =w 2 j-i) 

= pr(Ji = u 2 i-i\Xi = U2i-i) 

iy it pr(Ji = U2i-x) 

= pr(X i =u 2 i^i\Ji=u 2i -i) — — (5.4) 

pr(Xi = U2i-i) 

q t (l- a t -bj-i){(l - ai )(l - h)}- 1 
~ b i x 1 T. 

dibi 



(l-ai)(l-bi) 
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From equations (5.1)-(5.4), we obtain that 

pr(X l+ i = U2i\X\,. . .,Xi-i,Xi = u 2 i-i) 
1-ai-bi a t bi 

= 0i x — — — + 1 X 



(l-Oi)(l-6i) (l-Oi)(l-6i) 
b, 



1 - Oj 

Similar computations lead to 

rv iv v v \ (1 -6. t -Q t+ i)(l -a. t -&.,) 

(1 - Oj)(l - Oi) 

and 

, Y , v v v s g i+ i(l-a t -b t ) 

pr{X i+1 =u 2 i + 2\Xi,...,X i -i,X i = U2i-i) = -7- tt; rr> 

(1 - Oj)(l - h) 

which proves (4.3). □ 

Theorem 5.1 implies that ordered pivotal sampling shares the same second-order inclu- 
sion probabilities as Deville's systematic sampling. The computation of these probabilities 
is developed in Deville [4], and is reminded below. 

Theorem 5.2 (Deville [4]). Let k and I be two distinct units in U. If k and I are two 

non-cross-border units that belong to the same microstratum Ui, then 

TTfci = 0, 

if k and I are two non-cross-border units that belong to distinct microstrata Ui and Uj, 
respectively, where i < j , then 

TTfc; = 7r fe 7r;{l - c(i,j)}, 

if k = ki—i and I is a non- cross-border unit that belongs to the microstratum Uj where 
i < j , then 

TTfci =7T fc 7T;[l - 6j_x(l - 7T fe ){7T fc (l - l ) } ~ 1 c(l, j )] , 

if I = kj—i and k is a non-cross-border unit that belongs to the microstratum Ui where 
i < j , then 

TTfcz = 7r fe 7r ( {l - (1 - tt;)(1 - 6j_i)(7r i 6j_i)~ 1 c(i, j)}, 
if k = pi-i and I = Pj-i, where i<j, then 

TTfci = 7T k 7Tl[l - 6,-1 (1 - - 7Tfc)(l - ^{-Kk^lbj-l (1 - &i_ 1 ) } ~ 1 c(l, j )] , 

where c(i,j) = YViZi °i> c i = a ih{(l — — anc ^ C (M) = 1- 
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As noticed by Deville [4], it follows from Theorem 5.2 that many of the second-order 
inclusion probabilities are zero. As a result, no unbiased variance estimator may be found 
for the Horvitz-Thompson estimator. The search for variance estimators under reasonable 
model assumptions for the variable of interest y is a matter for further research. 

6. Interest of ordered pivotal sampling 

This is clear from Theorems 5.1 and 5.2 that ordered pivotal sampling induces a sampling 
design with a rather small entropy, since the second-order inclusion probabilities heavily 
depend on the order of the units in the population. If the maximization of entropy is a 
major concern, randomized pivotal sampling, where the list of the units in the population 
is randomly ordered before applying the pivotal method, should certainly be preferred. 
The main interest of ordered pivotal sampling lies in the gain of precision obtained from 
a stratification effect, if the ranking of the units in the population is well correlated 
to the variable of interest. In this sense, ordered pivotal sampling is similar in spirit to 
classical, ordered systematic sampling. However, systematic sampling can be particularly 
inefficient if the ordering is unappropriate, with regard to the variable of interest. Ordered 
pivotal sampling introduces more randomization in the sampling process, and should be 
more robust, in some sense, than systematic sampling. In the sequel, ordered pivotal 
sampling is compared to other sampling designs with respect to various criteria. 

To fix ideas, we consider the case of (i) equal inclusion probabilities iik = n/N, such 
that (ii) the population size TV is an integer multiple of the sample size n, and we note 
N = np. In this case, the microstrata Ui,i = 1, . . . ,n, are non overlapping with the same 
size Ni=p. We have 

Ui = {(i-l)p+l,...,(i-l)p + p}, (6.1) 

and ordered pivotal sampling amounts to stratified simple random sampling of size = 1 
inside each microstratum Ui. Also, it is well known that under the same assumptions (i) 
and (ii), systematic sampling amounts to simple random sampling of size to = 1 in the 
population G c = {gi, . . . , g p } of M = p clusters, where each cluster 

9j = {j,j+P,---,3 + (n-l)p} (6.2) 
contains Alj = n units. Let y denote some variable of interest, and let 

ty* = V) — ( 6 - 3 ) 

kes K 

denote the Horvitz-Thompson (HT) estimator of the total t y =Sj. g [/J/fc- 

Under conditions (i) and (ii), ordered systematic sampling and ordered pivotal sam- 
pling may be seen as particular cases of Markov chain designs (Breidt [2]). Let M be a 
doubly stochastic transition probability matrix, of size p. In a Markov chain design with 
matrix of transition M, a sample s = {R\,p + i?2, . . . , (n — l)p + R n } is selected, where 
i?i,...,i?„ is the Markov chain associated to M, with Ri being uniformly distributed 
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on {l,...,p}. Let I(p) denote the identity matrix of size p, and J(p) denote the square 
matrix of size p with all elements equal to 1. The use of the matrix of transition 

M p = p-M + (l- p )I(p) 7 

with p £ [0, 1] defines the category of compromise Markov chain designs (Brcidt [2]). The 
choice p = leads to ordered systematic sampling, while the choice p = 1 leads to ordered 
pivotal sampling. 



6.1. Entropies of sampling designs 

As a measure of randomness of a sampling design q{-), we use the entropy H{q) defined 
as 

#(<?) = -$>(*) log «(«), (6.4) 

sCU 

with OlogO = by convention. We have 

H(srs) = log A! - logn! - log(A - n)\ 

fc=0 v 7 

for simple random sampling, and 

H(sys) = log 

for ordered systematic sampling, see for example Tille and Haziza [14]. Some straight- 
forward algebra leads to 

H{ops) = ri log 

for ordered pivotal sampling. As a measure of comparison of entropy for two sampling 
designs g(-) and r(-), we may use the Kullback-Leibler divergence 

D(«||r) = 5>Wbsg4 

sCU ^ ' 

if the two sampling designs are such that r(s) = => q(s) = 0. We obtain 

D(sys\\srs) = ^ logf — J , 
fc=i ^ n ' 
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71 — 1 / _^ 



D(ops\\srs) = Y,los[y^), 



k=0 



(N 

D(sys\\ops) = (n — 1) log I — 

Both simple random sampling and ordered pivotal sampling clearly have much larger 
entropy than ordered systematic sampling. 

6.2. Maximum design-effect for sampling designs 

This is a standard fact that the variance of the HT-estimator under without-replacement 
simple random sampling is given by 

V srs {t yiT ) = N 2 ^-J-Sl (6.5) 

where / = n/N, S y = j^iJ^keuiVk ~ My) 2 and /j, y = jjY.keuVk- 0n the other nand , 
the variance of the HT-estimator under ordered pivotal sampling and assumptions (i) 
and (ii) may be written as 

V„ ps (U = JV 2 ^^^ ; , (6.6) 

i=l 

where = X) feet/ . iv* ~ Ms/*) 2 and N% = ~kl Y^keu, Vk- Finally, the variance of the 
HT-estimator under systematic sampling is then given by 

V sys {i y ^)=N 2 ^-^-Sl, (6.7) 
n n 



where 



2 P 
1 3=1 

with t y3 = £ feeG . Vk and m yj = hi In- 
ks a measure of risk of a strategy combining a sampling design q(-) and HT-estimation, 
we may use the maximum design-effect 

DMAX{q)=m^^^-, (6.8) 
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where C denotes the set of non-constant variables of interest (that is, containing all 
variables y such that Sy ^ 0). 

Theorem 6.1. Assume that conditions (i) and (ii) are satisfied. Then we have for or- 
dered pivotal sampling 

TV - 1 

DMAX(ops) = (6.9) 

TV — n 

and for ordered systematic sampling 

DMAX(sys) = n^—^-. (6.10) 

TV — n 

Proof. For any variable y, it follows from a standard analysis of variance that 

n _^ n 

S y = E 7v _ i ^ + E 7v - 1 ^' yi - 

i=l t=l 



so that 



E 5 y» - „ _ i S v 



p-1 y 

and the equality occurs if all the stratum means fi y i are equal. A joint application of 
(6.5) and (6.6) leads to 

Vapstfy*) TV-1 2/ 2= TV-1 

V srs {t yn ) -n(p-l) * TV-n' 

which gives (6.9). The use of an alternative analysis of variance leads to 

p , p 

n — 1 o x -> n , . 2 



s l = E ~fj~zrx a yj + E ^v - T (TOw ~ /iy)2 



where ct^. = ^ £ fceG . ~ m yi) 2 - This leads 



to 



2 n 2 TV-1 2 

Oy < : Jy, 

p — 1 n a 

and the equality occurs if the variable y is constant inside any cluster gj. By a joint 
application of (6.5) and (6.7), we have 

Vsysjiy.) n 2 N-l 2 2 N-l 

v srs (i v ^)- P -i m y' y /v-n' 

which gives (6.10). □ 
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If the sample size n remains small to moderate, equation (6.9) implies that DMAX 
tends to 1 in case of ordered pivotal sampling, if iV is sufficiently large. Even in the worst 
cases, ordered pivotal sampling will thus be competitive to simple random sampling. On 
the other hand, equation (6.10) implies that a strategy involving systematic sampling 
may be considerably more risky in some situations. 

6.3. Equality of treatment of variables for sampling designs 

As pointed out by a referee, the DMAX criterion considered previously is very stringent 
since referring to the worst possible variable for sampling designs. In case of ordered 
systematic sampling, this would be a cyclical variable whose period is equal to p = N/n: 
such a situation is usually unlikely to occur, except in particular situations. 

An alternative criterion studied by Devillc [3] and Qualitc [12] considers the equality of 
treatment of variables. For any sampling design q(-) with first-order inclusion probabilities 
Wk and second-order inclusion probabilities itklio) with k,l£U, let 

A(q) = [A kl (q)]k,ieu 

be its design variance-covariance matrix, with Aki{q) = ^ki(q) — ^k^i- Let < \\{q) < 
■ • • < Ajv(<?) denote the eigenvalues of A(q). For any variable y that lies on the unit sphere 
of the Euclidean norm (that is, such that ^2 kl£U yl — 1), we have 

Xl(q)<V p (i yv )<X N (q). 

Roughly speaking, the extreme eigenvalues give the extreme possible values for the vari- 
ance. Qualitc [12], page 50, proposed to measure the equality of treatment for variables 
in terms of minimization of the dispersion of the eigenvalues, denoted as 

N 

s( q ) = N-^{ x ^)-Mq)} 2 

fe=i 

with \{q) =N~ 1 Y l k=i X k(q)- Note that £f =1 A fc (g) =Tr(A(q)), where Tr(-) denotes 
the trace. For any sampling design q(-) with equal probabilities iik = n/N , this leads to 

\{q)=N- 1 Y J K k {l-Kk) 

keu 

_ P~l 

which will be simply denoted as A in the sequel. 

The ranking of the evaluated sampling designs with respect to this criterion is es- 
tablished in Theorem 6.2. Clearly, ordered pivotal sampling tends to treat the variables 
more equally than ordered systematic sampling. To demonstrate Theorem 6.2, we need 
the following lemma. 
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Lemma 6.1. Assume that the sampling design q(-) is performed with equal probabilities 
7Tfc =n/N, and that conditions (i) and (ii) are satisfied. Assume that A(q) has only two 
eigenvalues and \+{q) > 0, with multiplicities No(q) and N — No(q) , respectively. Then: 



S(q) = 



No(q) _- x2 _ 



N-N (q) 



(6.11) 



Proof. Since q(-) has only one strictly positive eigenvalue A+(g) with multiplicity N ■ 
N (q), we have X + (q) = N _^ o{q ) X. Then 



JV 



N5(q)=Y,{Ml)~W 



fe=l 



= N {q)X 2 + {N - N (q)} 



N 



N-N (q) 



A- A 



N- 



No(q) 



-A . 



N-No(q)"' □ 
Theorem 6.2. Assume that conditions (i) and (ii) are satisfied. Then we have 

5(srs) < S(ops) < S(sys). (6.12) 

Proof. It may be easily shown (see, e.g., Devillc [3], page 120) that in case of simple ran- 
dom sampling, A(srs) has only two eigenvalues and X+(srs) = , with No(srs) = 1. 

Let I(n) denotes the identity matrix of size n, and J(n) denotes the square matrix of 
size n with all elements equal to 1. After some algebra, we obtain from the definition of 
ordered pivotal sampling that 



A (ops) = 



/Ax 
Ai 

V ••• 







\ 



'■. 

o aJ 



J(n) ® Ai 



where Ai =p~ 1 {I(p) — p~ 1 J(p)} and ® denotes the Kronecker product. It follows that 
the N eigenvalues of A(ops) are given by the products of the eigenvalues of I(n) and Ai 
(see for example Theorem 1 in Magnus and Neudecker [10], page 28). The eigenvalues of 
Ai are p~ x and with multiplicities p—1 and 1, respectively. Consequently, A(ops) has 
only two eigenvalues and A+(ops) =p _1 , with No(ops) — n. 

Similarly, we obtain from the definition of ordered systematic sampling that 



/A, 



AA 



A(sys) 



J(n) <g> Ai 



\Ai 



Ar/ 
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Table 2. Values of three variables of interest in the generated population 



Unit 


l 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


12 


yi 


10 


10 


10 


15 


45 


45 


50 


50 


60 


60 


60 


65 


2/2 


15 


45 


10 


60 


60 


50 


45 


65 


10 


50 


10 


60 




10 


45 


60 


15 


50 


65 


10 


50 


60 


10 


45 


60 



Since the eigenvalues of J(n) are and n with multiplicities n — 1 and 1, respectively, 
A(sys) has only two eigenvalues and X + (sys) = np^ 1 , with N (sys) = N — p + 1. 

Equation (6.11) implies that S(q) increases as N (q) increases. Clearly, N (srs) < 
No(ops), and from the identity (N — p + 1) — n = (p — l)(n — 1) > 0, we obtain that 
N (ops) < N (sys) so that the result follows. □ 



6.4. Some numerical results on a small population 

To investigate on the properties of considered sampling algorithms, we considered a 
small example. We first generated a finite population of size N = 12, containing three 
variables of interest, y±, y^ and y%. Table 2 shows the values for the three variables of 
interest. The variable y\ is highly correlated to the order of the units in the population, 
on the contrary to variable y2- The variable 1/3 exhibits a particularly unfavorable case 
for systematic sampling. 

We considered equal probability sampling of size n = 2 (respectively, n = 4) by means 
of six sampling designs: simple random sampling without replacement (SRS), ordered 
systematic sampling (SYS), compromise Markov chain design with p — 0.25 (CMC25), 
p = 0.50 (CMC50), p = 0.75 (CMC75), and ordered pivotal sampling (OPS). As a measure 
of variability of the HT-estimator t y7r for a sampling design q(-), we considered the design- 
effect (DEFF) given by 

DEFF = Vq{ - ty ^ , (6.13) 
V (t V 

v srs x^yir ) 

where the variances are computed by means of formulas (6.5)-(6.7) for SRS, OPS and 
SYS, and from formula (1) in Breidt [2], page 66, for compromise Markov chain designs. 
Table 3 shows DEFF for the five strategies. As could be expected, the CMC25, CMC50 
and CMC75 give compromise results between SYS and OPS. Also, it is clear from Table 3 
that both OPS and SYS lead to a subsequent reduction of variance for variable j/i, with 
DEFF ranging from 0.17 to 0.50 and OPS performing significantly better. The OPS 
strategy is essentially similar to SRS for the variable which is poorly correlated to the 
order of the units in the population, while SYS may be much worse (DEFF = 1.39 for 
n = 2) or much better (DEFF = 0.36 for n = 4). Finally, we obtain for the variable 7/3 a 
considerable loss for SYS, while the loss is more limited for OPS with DEFF = 1.10 for 
n = 2 and DEFF = 1.36 for n = 4. 
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Table 3. Design-effect for three variables of interest and five strategies in the generated popu- 
lation 





Sample size 


n = 2 




Sample size n = 4 






2/i 


V2 


2/3 


2/i 


V2 


2/3 


SYS 


0.50 


1.39 


2.18 


0.27 


0.36 


5.44 


CMC25 


0.46 


1.31 


1.91 


0.24 


0.61 


3.94 


CMC50 


0.43 


1.24 


1.64 


0.21 


0.76 


2.81 


CMC75 


0.39 


1.17 


1.37 


0.19 


0.85 


1.97 


OPS 


0.35 


1.10 


1.10 


0.17 


0.95 


1.36 
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